Chaos of the Relativistic Parametrically Forced van der Pol Oscillator 
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A manifestly relativistically covariant form of the van der Pol oscillator in 1+1 dimensions is studied. 
We show that the driven relativistic equations, for which x and t are coupled, relax very quickly to 
a pair of identical decoupled equations, due to a rapid vanishing of the "angular momentum" (the 
boost in 1 + 1 dimensions) . A similar effect occurs in the damped driven covariant Duffing oscillator 
previously treated. This effect is an example of entrainment, or synchronization (phase locking), of 
coupled chaotic systems. The Lyapunov exponents are calculated using the very efficient method of 
Habib and Ryne. We show a Poincare map that demonstrates this effect and maintains remarkable 
stability in spite of the inevitable accumulation of computer error in the chaotic region. For our 
choice of parameters, the positive Lyapunov exponent is about 0.242 almost independently of the 
integration method. 

PACS numbers: 05.45. +b, 04.20.Cv 
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The chaotic behavior of a covariant relativistic gener- 
alization of the classical Duffing oscillator in 1 + 1 di- 
mensions has been studied recently [0J|] . This system is 
completely integrable, since both the total "mass" (the 
value of the invariant generator of dynamical evolution) 
and the hyperbolic "angular momentum" (the generator 
of Lorentz transformations in 1 + 1 dimensions) are con- 
served. Under non-stationary perturbation, including a 
friction term as well, however, the stable and unstable 
orbits in the phase space separate and cross (infinitely 
many times). ^] 

Some general features of the resulting chaotic orbits 
were discussed in ref. [Q , such as the phenomena of pe- 
riod doubling and the apparent existence of a strange 
attractor ||. In ref. analytic solutions were stud- 
ied for the unperturbed problem in the neighborhood of 
the separatrix, and it was shown that the perturbative 
Mclnikov criterion Q| for homoclinic chaos was satisfied. 
It was found that along with chaotic behavior in space, 
there is chaos in time in such a system, with perhaps 
profound implications for our perception of the nature of 
dynamical processes. 

In this paper, we study a manifestly covariant dynam- 
ical system with no Hamiltonian structure, i.e., a covari- 
ant form of the van der Pol oscillator. The model was 
originally introduced by van der Pol || to describe the 
triode electronic oscillator, and used by van der Pol and 
van der Mark as a model for heart rhythms ||. The van 
der Pol equation was introduced before chaos was sys- 
tematically defined or studied, and it became one of the 
first examples to be analyzed for chaotic behavior and 



*One can think of the friction term as the result of radiation 
reaction for a charged oscillator, i.e., the damping of a radi- 
ating dipole which depends only on the relative motion of the 
charges 



associated entrainment phenomena j7J-10|. This system 
carries intrinsic (coordinate dependent) friction; if this 
term vanishes, one obtains the structure of an ordinary 
harmonic oscillator. The method of Melnikov is not ap- 
plicable for the van der Pol system, since there is no ho- 
moclinic orbit associated with the unperturbed problem. 
We shall therefore calculate the Lyapunov exponents to 
demonstrate the chaotic behavior [Q. 

It was observed in the case of the relativistic Duffing 
problem that the x and t modes, viewed as two coupled 
systems, could converge, in the presence of friction, to a 
single Duffing type oscillator (the "angular momentum" 
converges to zero) . This convergence is an example of the 
entrainment (sometimes called "phase locking") of two 
coupled chaotic systems. We observe a similar behavior 
in the relativistic van der Pol equations. It is conceivable 
that nonconservative relativistic electronic systems have 
similar properties, providing examples of fundamental 
interest for the entrainment, or synchronization (phase 
locking) properties of such chaotic systems . 

As for the Duffing oscillator, we add a driving force 
that conserves the "angular momentum" , so that there 
is minimal distortion introduced by such a term. We find, 
as for the Duffing oscillator, that the time and space 
modes coalesce, with resulting loss of independent de- 
grees of freedom. The resulting one-dimensional system 
does not coincide precisely with the usual nonrelativistic 
case; the forcing term appears proportional to the de- 
pendent variable (x(r) or t(r)). This term must be then 
shifted to the non-linear coefficient of the dependent vari- 
able in the van der Pol equation (previous studies have 
considered driving terms in the coefficient of the x jl2| ) . 
It was therefore necessary to investigate the chaotic struc- 
ture of the resulting system independently of previous 
studies. 

Although the chaotic behavior of the usual van der Pol 
system is difficult to observe, we found that in our case, 
the usual methods utilizing, for example, Poincare maps, 
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clearly demonstrate the more complex chaotic character 
of this system. Moreover, we employ a method given 
by Habib and Ryne |ll| for the study of the Lyapunov 
exponents to confirm the chaotic nature of the system; 
this method is very effective in studying a system of this 
type, i.e., a system for which the local linearization can 
be described in terms of a Hamiltonian form. 

The rapid relaxation of the dissipative relativistic sys- 
tem to a system closely related to the nonrelativistic form 
is of interest in itself in providing an example of how dissi- 
pation can induce an effectively nonrelativistic dynamics 
(through cntrainmcnt of the time and space modes). A 
similar phenomena occurs, as mentioned above, in the 
dissipative form of the Duffing oscillator under certain 
conditions (!]]. 

Following the classical relativistic mechanics of Stueck- 
elberg |Q , and its extension Jl5| to the many body prob- 
lem, one can define a relativistic invariant evolution func- 
tion K (analogous to the nonrelativistic Hamiltonian), 
which generates "Hamilton-like" equations: 
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dK 
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where and are the energy-momentum and space- 
time coordinates (i = 1 ... N, fj, = 0,1, with metric -, +) 
and t is a universal invariant parameter for motion in 
27V-dimensional phase space. 

Consider the two body problem with a potential that 
depends on the invariant distance between them, for 
which the corresponding evolution function, K, may be 
written as: 



K = 



V(p 2 ), 



2Mi 2M 2 
where p 2 is the Poincare invariant 

P 2 = (^i ~ x&ixin ~ x %ll ) 
= {xi - x 2 f - fa - t 2 f 



(2) 



(3) 



It is possible to write the two body problem as a function 
of center of mass and the relative motion between the 
bodies (in a way similar to nonrelativistic mechanics): 
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where 



K = K, 
K — 

K re l = 

M = Mi + M 2 



Pi + P2 



K re l 



(4) 
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FIG. 1. The angular momentum, Af 01 (r), of the typical 
cases where there is a) exponential growth followed by an os- 
cillatory exponential decay of M 01 (r) where M$ 1 Mo 1 > 0, b) 
oscillatory exponential decay of M 01 (r) where Mq 1 Mq 1 < 



The equations of motion in the new variables are 

dK rel 



dx M dK re i 
dr dp^ 

dX^ dK c . m . 



dp^ 
~dV 



(5) 



dP^ 



dr dP^ dr 



= 0; 



the system separates to two equations describing the cen- 
ter of mass and the relative motion of the system. 
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FIG. 2. The (y, z) plane of fig. [j]a and b. The inset in fig. 
^> shows an enlargement of the turning point from the initial 
point to the final linear curve. 

The equation for the nonrelativistic externally forced 
van dcr Pol oscillator is §,§ [] : 

x + a(x 2 — l)x + kx — g cos tot, (6) 

where the right hand side corresponds to an external 



t Jackson M has shown that a self-stabilizing radiating sys- 
tem can be described by this equation, where x corresponds 
to dx jdt of the original system. The quantity dx^' jdr oc- 
curs in the relativistic form of such a system, where dr is an 
invariant. The same substitution would then be valid in this 
case. 



driving force. Although there is no Hamiltonian which 
generates this equation, we may nevertheless consider its 
relativistic generalization in terms of a relative motion 
problem, for which there is no evolution function K. 

In this paper, we study the relativistic generalization 
of Eq. (0) in the form 



+ a{p 2 - Vjx^ + x^ = 0, 



or, in terms of components, 



x + a(x 2 — t 2 — l)x + kx = 
i+a(x 2 - t 2 - l)i + kt = 0, 



(7) 



(8) 
(9) 



where x, t are the relative coordinates of the two body 
system. The term x can be understood as representa- 
tive of friction due to radiation, as for damping due to 
dipole radiation of two charged particles, proportional to 
the relative velocity i M in the system j|. To study the 
existence of chaotic behavior on this system, we add a 
driving forcej[] to the system (it already contains dissipa- 
tion intrinsically) in such a way that it does not provide 
a mechanism in addition to the dissipative terms for the 
change of "angular momentum" (this quantity is con- 
served by eq. ([]) with a = 0) 



M' 



01 



m{tx — xt). 



(10) 



To achieve this, we take a driving force proportional to 
x^, so that (||) becomes 

x + a(x 2 — t 2 — l)x + kx = fxcosujT (11) 
i + a(x 2 -t 2 -l)i + kt = ft cos lut. (12) 

The derivative of the angular momentum can be evalu- 
ated in terms of these equations as 



dM° 



-a(x 



1)M 



01 



(13) 



In contrast to exponential decay of the angular momen- 
tum that was shown in Duffing oscillator W , one cannot 
see clearly whether eq. (|13| ) reflects such decay or not. 
Exponential decay of M°^neans that after a relatively 
short time M 01 goes to zero, and that x becomes propor- 
tional to t (x = (3t). Thus, in this case, the two coupled 
equations would actually reduce to one equation. In or- 
der to simplify the investigation of eqs. ( Ill ) and dl3), 
one can introduce two symmetric equations by subtract- 
ing and adding them, to obtain 



y + a(yz - l)y + ky = fy cos ujt, 
z + ot(yz — \)z + kz = fzcosujT, 



(14) 
(15) 



'We understand 'if 1 to be proportional to force, from its form 
in the framework of the Hamiltonian theories. 
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where y — x — t and z = x + t. Eqs. (|14|), ( |15| ) are 
completely symmetric and it is clear that if one starts 
from initial conditions {yo,yo) = (3(zo,zo) the equations 
behave identically (up to a multiplicative constant). In 
this case, the angular momentum, M 01 = ^(zy — yz), is 
zero, and hence y — /3z. The derivative of M in 
new variables is, 



dM 



01 



dr 



-a(yz - 1)M° 



(16) 



For the system (|T|), (|l|), M 01 appears to converge, in 
actual computations, to zero independently of initial con- 
ditions. In fact, there are three typical behaviors of M 01 
before it reaches zero: (a) An exponential growth (or os- 
cillatory exponential growth) of |M 01 | followed by expo- 
nential decay (or oscillatory exponential decay) to zero, a 
scenario which occurs when M^M® 1 > 0. (b) An expo- 
nential decay (or oscillatory exponential decay) of |M 01 | 
to zero, which occurs when M 01 ^ and M^Mg 1 < 0. 
(c) M 01 is zero initially and remains zero. 

In fig. |l| we present the behaviors of cases (a) and 
(b). The numerical calculation of eqs. (|l4j), ([l5|), were 
performed using the adaptive time step fifth-order Cash- 
Karp Runge-Kutta method |T(|, and by fixed time step 
forth-order Runge-Kutta method jl6| . We choose m = 1 
for the mass and to = 2.466, a = 1, k = 1 and / = 10. 
The initial conditions used in fig. [tja were yo = 1, and 
yo = ZQ = z = -1 (M 01 = 1 and Af 01 = 2), and 
it is easy to observe the initial exponential growth and 
the succeeding oscillatory exponential decay to zero. In 
fig. [j]b a typical behavior of case (b) is shown. The 
initial conditions are yo — 3, yo — 2, zo — | and io = 
(Mq 1 = i and A/g 1 — j), and the typical behavior of an 
exponential oscillatory decay to zero is shown. If fig. |2|a 
and ||b we show the (y, z) plane of fig. [j]a and |l|b which 
the solutions converge to a linear curve after a short time 
and then continue with oscillations on this line. Hence 
we also see in the phenomenon of entrainment in the van 
der Pol case. 

After the angular momentum reaches zero there is no 
reason to investigate the two coupled eqs. (|l4j), (15), 



which actually become uncoupled since the variables are 
proportional. Thus, one is left with the equation (we 
henceforth replace r by the familiar designation t appro- 
priate to the nonrelativistic case): 



+ a(x 2 — l)x + kx = fx cos(<x>i). 



(17) 



A class of equations including cq. jl7|) were studied 
by Schmidt and Tondl [jlTj for the case for which both 
damping and nonlinear terms, as well as the homoge- 
neous driving term are taken as small (as can be seen by 
scaling the dependent variable); they used perturbative 
techniques to study the instability. Cicogna & Fronzoni 
jl8| have studied the modification of chaos of the Duffing 
oscillator using an additional forcing term which can be 
proportional to similar structure occurs in ref. [pi, 



as mentioned above. Such a forcing term occurs in the 
Mathieu equation (as a time-dependent frequency) [fT9| . 
Our study is not restricted to small coefficients; it utilizes 
numerical methods of wider applicability. 

In the following we will show the chaotic properties of 
eq. ( |l7j ) , including the complex shape of the trajectory 
itself (which is reflected by a continuous frequency spec- 
trum), and the fractal form of the Poincare map, which in 
our case is just (x(t n ),x(t n )) where t n = t a + nT (T = ^j- 
and n is an integer number larger then 0). The chaotic 
behavior is verified by the existence of a positive Lya- 
punov exponent. 

There exist several methods for calculating Lyapunov 
exponents from a set of first order differential equations 
po| , pT| . However, we preferred a new method developed 
by Habib and Ryne Jl5| ], which is based on a technique 
using symplectic matrices. The symplectic computation 
of Lyapunov exponents is applicable whenever the lin- 
earized dynamics is Hamiltonian (which true for our sys- 
tem). The basic advantage of the symplectic method is 
that it avoids the renormalization and reorthogonaliza- 
tion procedures necessary in the usual techniques, and 
thus, leads to fast and accurate results (the method 
avoids an additional error which can be caused by the 
renormalization and reorthogonalization procedures). 

According to the method of Habib and Ryne one lin- 
earizes the system, as a first step in the computation of 
Lyapunov exponents. In our case, the linearization of eq. 
©is 

<5 + a{xl - 1)8 + {2ax xo +k- f cos ut)S = 0, (18) 

where S — x — xo and xo is the fiducial trajectory 
Following Habib and Ryne, we change the variable to 
A = <5cxp(— g(t)), where g(t) = —§{xl — 1), so that eq. 
( |l8| ) becomes, 

A - (g(t) + g(t) 2 - k + f cos ut)A = 0. (19) 

In order to find the Lyapunov exponents it is necessary 
to solve two first order differential equations, |l3| 

dfi 1 

^ = -( S22 - S n)cosa 

da . . 

— = sn + s 2 2 - (S22 - six) smacoth/i. (20) 

dt 

The Lyapunov exponents are then ±/i/i, where t is suf- 
ficiently large. These equations are obtained as follows. 
In the general case of a Hamiltonian of quadratic form, 



^ 2m 

H(Z,t) = ~ SijZtZj, 



(21) 



<-j=i 



where Z = (qi, q 2 , . . . , q m ,Pi,P2, ■ ■ ■ ,Pm), the canonically 
conjugate coordinates and momenta, the equation of mo- 
tion for Z(t) = M(*)Z(0) implies that 



—MM = JSMM - MMSJ, 
dt 



(22) 
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where 



and S 



J 



1 

-1 



The equations (20) are then obtained for 



space is represented in four dimensions the mapping of 
chaotic regions are quite difficult and it is necessary to 
(23) focus on one or two parameters; we choose a to be the 
varying parameter since it controls the nonlinearity of 
the system. 



a system with one degree of freedom (in the case that 
there is no term of the type qp in the Hamiltonian) , by 
recognizing that the most general form of the symplectic 
matrix M is given by 



M 



g/i(B2 cos a+B 3 sin a) ^>B\ 



(24) 



where B\ — ia\, B2 — 02 and B3 — (73, and {(Ji\ are the 
Pauli matrices. 

In the case of eq. (p5[), the are : 



six 



-{g{t) + aitf - k + fcostot) 



= (XXqXo 

S22 = !• 



2 



l) 2 + k — / cosuit, 



(25) 



After finding the Lyapunov ex- 

ponents, X± = ± lim^oo [i(t)/t, of eq. (19) using eqs. 
( p0| ) , one must find the Lyapunov exponents of the origi- 
nal system, eq. ( |l8| ) . By recognizing that asymptotically 
eq. ( |l9| ) is A± = 5± exp(— g(t)) = exp X±t, one can then 
find the Lyapunov exponents of eq. (Il8]) fl3|] : 



X ±= Hm Ug{t) ± fi(t)). 

t^oc t 



(26) 



It is possible to summarize the Lyapunov exponents 
calculation of our system (eq. (|l7|)) by five first ordered 
differential equations: 



<1xq 
~~dt 



Po, 



= -a(xl - l)p - kxo + fx cos ut, 

dg _ at, 2 
dt~ 2 (X ° 



0.5 



0.0 



-0.5 



-1.0 




0.0 
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10.0 

a 



15.0 



20.0 



1 



(27) 
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1 



(1 — ax a p a + — (xq — l) 2 — k + f cos ut) cos a, 



-50 



(b.) 
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da 

— = 1 + ax p 
dt 



— (xq - l) 2 + k - f cos u)t 



10.0 

a 



15.0 



20.0 



a 2 2 2 FIG. 3. Lyapunov exponents as a function of a for a) 

(1 — axopo H — (xq — 1) — k + f cos ut) sin a coth /x, positive Lyapunov exponent, b) negative Lyapunov exponent. 



and the Lyapunov exponents can be found by eq. (26). 
The first two equations calculate the fiducial trajectory. 
The third relation is needed for the calculation of the 
Lyapunov exponents. The last two equations are just 
the explicit form of eqs. (p0|). 

In order to locate some parameter values which lead 
to chaotic behavior, one must map the positive Lya- 
punov exponent in parameter space. Since our parameter 



In fig. H we show Lyapunov exponents as a function 
of a (the other parameter values are lo — 2.466, k = 1, 
/ = 5). We draw the Lyapunov exponents after 200000 



§The parameter space is actually three dimensional since 
one can avoid the parameter k in eq. ([Tt]) by introducing a 
rescaled time s = \/\k\t. 
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time units. The larger Lyapunov exponent is shown in 
fig. 0a, and it is easy to see that it becomes positive 
in a small region of the range of a. One can also iden- 
tify some periodicity in x+( a ) multiplied by a decaying 
envelope. This "periodicity" gives a clue about some 
intrinsic frequency which resonates with the "external" 
frequency, u>. However, it seems that the gap between 
the chaotic regions is not exactly constant. A negative 
Lyapunov exponent is shown in fig. |^b. As expected, 
the absolute value of the negative Lyapunov exponent is 
much larger than the positive one. The negative Lya- 
punov exponent decreases approximately linearly with a 
(not including some regions parallel to chaotic regions 
with positive Lyapunov exponent) since a controls the 
dissipation of the system. 

One of the fundamental questions in the numerical so- 
lution of chaotic systems is the validity of the numerical 
results. Since in every numerical integration there is an 
accumulative error, the integrated trajectory separates 
from the real trajectory in a very short distance. Accord- 
ing to chaos theory, a nearby trajectory can generate a 
completely different trajectory (nearby trajectories actu- 
ally separate exponentially), thus, one cannot guarantee 
the accuracy of the numerical calculation. The numerical 
trajectory is method dependent, and more than this, is 
accuracy dependent (if one changes the integration ac- 
curacy the trajectory behaves differently). One should 
expect, moreover, sensitivity to the integration method 
(or, equivalently, sensitivity to accuracy of the integra- 
tion) if the system is chaotic, a fact which can used as 
an additional sign of chaotic behavior. Thus, it is clear 
that one cannot rely on the numerical results. The ques- 
tion is whether one can rely on the global characteristic 
values, such as, Lyapunov exponents. It seems that the 
answer is that one can indeed trust the global values, al- 
though the local behavior is not assured Q It seems that 
although the trajectories behave differently they appear 
to "move" on the attractor; they have the same global 
characteristic exponents. The influence of this sensitivity 
can be seen in fig. ||a. One recognizes easily that there 
is a clear smooth behavior of x+( a ) when it is less then 
zero (the system is not chaotic and nearby trajectories do 
not separate exponentially); once the exponent becomes 
positive x+( a ) becomes less smooth. 




0.0 50000.0 100000.0 150000.0 

t 




0.0 50000.0 100000.0 150000.0 

t 

FIG. 4. Convergence of Lyapunov exponents to a constant 
value, as a function of t, using four different integration meth- 
ods for a) positive Lyapunov exponent, b) negative Lyapunov 
exponent. 



** Greene |22j pointed out that there are huge numerical er- 
rors in calculating orbits in stochastic regions, yet these errors 
do not seem to change the fundamental character of the orbit. 
In his appendix he shows analytically that the error propa- 
gates anisotropically according to the direction of maximal 
instability. The stability transverse to these region was aston- 
ishingly well preserved and thus the boundaries of the region 
were visually stable regardless of instability along the orbit. 
Our results are in accordance with Greene's observations. 
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FIG. 5. 
oscillator. 



The attractor of relaxed relativistic van der Pol 



In order to verify the claim that the Lyapunov expo- 
nents are the same for different integration techniques, 
we compare the results between four different methods: 
1. Bulirsch-Stoer method using Richardson extrapola- 
tion jl6| (denoted by bsstep); 2. adaptive stepsize fifth 
Runge-Kutta method using Cash-Karp parameters 
(denoted by rk5); 3. fourth Runge-Kutta method 
(denoted by rk4); and 4. Adams integration method 
We choose a — 1.149 since it gives the largest positive 
Lyapunov exponent as seen in fig. The Lyapunov ex- 
ponents are shown in fig. ^ and it is clear that all methods 
lead to approximately the same exponents; \+ = ~ 0.242, 
X- =~ -4.286. 

One can find other signs of chaotic behavior with the 
parameters used in fig. || (u — 2.466, k = 1 and / = 5), 
and taking a = 1.149. In fig. |^ the real phase space, 
(cj), x,p), where (f> — tot mod 2ir, is shown. It can be seen 
clearly that the trajectory is very complicated and does 
not return to itself. The phase space lies on a surface with 
some width; a fact which suggests a noninteger fractal 
dimension. Using the Kaplan and Yorke conjecture [£4| , 
the fractal dimension would be : 



Dl = J 



IVnl 



0.242 
4.286 



2.06 



(28) 



where j is defined by the condition that ^i=i K > 
and Yli=i < (the Lyapunov exponent in the time 
direction is, of course, zero). 



FIG. 6. The complex behavior of x(t) (we have used the 
notation t instead of r to describe the decoupled system). 

No clear periodicity is seen in x(t) , as shown in fig. |^. 
The Fourier transform of x(t) (fig. |7J) shows a continuous 
frequency spectrum, which is a sign of chaotic behavior. 
In addition, there is an approximate linear decay of the 
logarithm of the power spectrum, a fact which is con- 
nected to the exponential decay of the autocorrelation 
function [^5| (we also checked the autocorrelation itself 
and its envelope indeed decays exponentially). One sees 
a peak in the neighborhood of 1, which is the frequency of 
the undamped and unforced oscillator. One would expect 
resonant behavior in the neighborhood of half the driving 
frequency (O = on the basis of a comparison with the 
Mathieu equation fll9| . In our case, there is a apparent 
resonant region around this frequency £1 ~ 1.233. 

Despite the difficulties mentioned above, we succeeded 
to draw a Poincare map for eq. (]T7]) . The map in this case 
is just the points on the plane <fi = 0, and can be written 
as (x(t n ),p(t n )), where t n — ^-n and n is a nonnegative 
integer. The map is drawn using two different methods 
: rk4, and the Adams method; one obtains the same 
mapping structure. This is a surprising result since, as 
mentioned before, one can not be sure of the accuracy 
of the trajectory x(t) due to accumulation error which 
leads to a different trajectory behavior. The map shown 
in fig. U shows smooth curves with stretching and folding 
phenomenon (more fine structure is seen in the inset). 
The map is antisymmetric, and obviously has a fractal 
dimension w 1.06 (following eq. ( p8| ) and ignoring the 
time direction). 
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O.O 1 5.0 10.0 15.0 20.0 

Q. 

FIG. 7. The power spectrum of fig. []. The dashed line 
indicates the free oscillator frequency, i.e., Q = yk = 1. 

The cntrainmcnt phenomenon observed between space 
and time chaotic modes, both for the Duffing 0] and 
van der Pol case, for which x and t appear to approach 
exponentially to a smooth linear relation, even though 
both x and t remain chaotic, presents a mechanism for 
the rapid convergence of damped relativistic systems to 
their classical limit. Underlying this smooth behavior of 
x(t), one might expect to see a characteristic property of 
the radiation field reflecting the chaotic behavior of both 
x and t, in such systems. 

We further remark that the stability of the Poincare 
plots of fig. ^ indicates that computational deviations of 
x are compensated to reach a smooth curve by deviations 
in t. This stability is analogous to the stability of the 
calculation of the Lyapunov coefficients along a fiducial 
curve that is calculated with an inevitable error. Both 
the Poincare plots and Lyapunov calculations appear to 
reflect properties of the attractor rather that the local 
computed orbit. 
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FIG. 8. The Poincare map of the relaxed relativistic van 
der Pol equation. The inset shows an enlargement of the map. 
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